from pulp import *
import pandas as pd
import matplotlib.pyplot as plt 

import warnings
warnings.filterwarnings("ignore")

import pandas as pd
hourly_demand = pd.read_csv('./demands.csv', parse_dates=['hora'], index_col='hora').sort_index()
hourly_demand

hourly_demand.plot(figsize=(16, 4), legend=False, ylabel='demanda (MW)', xlabel='día y hora');

daily_demand = hourly_demand.iloc[:24]   # los 24 primeros registros
daily_demand.T   # trasponer para que no ocupe tanto en pantalla

daily_demand.plot(marker='o', legend=False, figsize=(8, 4), ylabel='demanda (MW)');

plants = pd.read_csv('plants.csv', index_col='central').sort_index()
plants

plant_costs = pd.read_csv('./plant_costs.csv', index_col='combustible').sort_index()
plant_costs

health_costs = pd.read_csv('./health_costs.csv', parse_dates=['hora'], index_col=['hora', 'central']).sort_index()
health_costs

health_costs.reset_index().pivot(index='hora', columns='central', values='coste_salud').plot(figsize=(16, 4), ylabel='coste sanitario', xlabel='día y hora');

daily_health_costs = health_costs.iloc[:72]   # 3 centrales, cada una de ellas 24 horas -> 72 filas
daily_health_costs.reset_index().pivot(index='hora', columns='central', values='coste_salud').plot(marker='o', figsize=(8, 4), ylabel='coste sanitario');

# Definición del problema
prob = LpProblem("Optimizacion_Energia", LpMinimize)

C = list(plants.reset_index()['central']) # Conjunto de centrales C
H = list(daily_demand.index) # Conjunto de horas que componen un día H

def crear_variables_decision(C, H):
    # Variables de decisión
    e_ch = LpVariable.dicts('MWh', (C, H), lowBound=0, cat=LpContinuous) # Energía producida por la central c durante la hora h
    x_ch = LpVariable.dicts('X', (C, H), cat='Binary') # Variable binaria que indica si la central c esta encendida o no en la hora h
    gamma_ch = LpVariable.dicts('gamma', (C, H), cat='Binary') # Variable binaria que indica si la central c está encendida (1) o no (0) en la hora h.
    return e_ch, x_ch, gamma_ch

e_ch, x_ch, gamma_ch = crear_variables_decision(C, H)

alpha_values = [0.8 if c in ['D', 'E'] else 0.01 for c in C]
df_alpha = pd.DataFrame({
    'central': C,
    'alpha': alpha_values
}).set_index('central')

alpha = df_alpha.squeeze() #Porcentaje de producción de la central c


beta_values = [0.2 if c in ['D', 'E'] else 0.25 if c in ['A','B', 'C', 'J'] else 1 for c in C]
df_beta = pd.DataFrame({
    'central': C,
    'beta': beta_values
    }).set_index('central')

beta = df_beta.squeeze() #Porcentaje de cambio mínimo de producción de energía por central

#Producción máxima
E_c = plants.reset_index()[['central','producción']].set_index('central') #Energía máxima de energía para una central c

## Costo del combustible
P_c = plants.join(plant_costs, on= 'combustible').reset_index()[['central','coste_comb']].set_index('central') # Unimos las tablas y seleccionamos las columnas necesarias
P_c.index.name = None # Unimos las tablas y seleccionamos las columnas necesarias
for h in H:
    P_c[h] = P_c['coste_comb']
P_c = P_c.drop(['coste_comb'], axis=1) # Costo de combustible de la central $c$.

## Costos operativos
O_c = plants.join(plant_costs, on= 'combustible').reset_index()[['central','coste_oper']].set_index('central')
O_c.index.name = None # Unimos las tablas y seleccionamos las columnas necesarias
for h in H:
    O_c[h] = O_c['coste_oper']
O_c = O_c.drop(['coste_oper'], axis=1) # Costo de operativo de la central $c$ por hora.

## Costos sanitarios 
S_ch = health_costs.reset_index().sort_index() # Costo sanitario por hora de la central $c$ en la hora $h$.
S_ch = S_ch.pivot(index='hora', columns='central', values ='coste_salud') # Pivotamos el DataFrame 'S_ch' para reorganizarlo y facilitar su análisis.

non_carbon_plants = ['D', 'E', 'F', 'G', 'H', 'I', 'J']

# El método reindex se utiliza para agregar las plantas sin costos de salud al DataFrame.
# El parámetro 'fill_value' se establece en 0, lo que reemplaza todos los valores faltantes con 0.
S_ch= S_ch.reindex(columns=S_ch.columns.tolist() + non_carbon_plants, fill_value=0)

D_h = daily_demand # Demanda total de energía en la hora $h$.

ON_c = plants.join(plant_costs, on= 'combustible').reset_index()[['central','coste_on_off']].set_index('central')
ON_c.index.name = None # Unimos las tablas y seleccionamos las columnas necesarias
for n in H:
    ON_c[n] = ON_c['coste_on_off']
ON_c = ON_c.drop(['coste_on_off'], axis=1) # Costo de encender/apagar de la central $c$

def crear_funcion_obj(prob, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, C, H):
    # Función: Costo de combustible
    comb = (lpSum([P_c[h][c] * e_ch[c][h] for c in C for h in H]))

#Función : Costo operacional
    op =  (lpSum([O_c[h][c] * x_ch[c][h] for c in C for h in H]))

# Función: Costo sanitario
    san = lpSum([S_ch[c][h] * e_ch[c][h] for c in C for h in H])

# Función costo de encender o apagar una central (Se define para su uso a posteriori, puesto que para este ejercicio nunca se aplicará)

    on_off = (lpSum([ON_c[h][c] * gamma_ch[c][h] for c in C for h in H]))

# Definiendo la función objetivo como:
    
    prob+= comb + op + san + on_off
    
    return prob

prob= crear_funcion_obj(prob, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, C, H)

def anadir_restriccion_demanda(prob, e_ch, D_h, C, H):
    ## Satisfacción de la demanda
    for h in H:
        prob += (lpSum(e_ch[c][h] for c in C) >= D_h.loc[h][0]) # Satisfacción de la demanda

    return prob

prob= anadir_restriccion_demanda(prob, e_ch, D_h, C, H)

def anadir_restriccion_xch(prob, x_ch, C, H):
    # Todas las centrales están encendidas en todas las horas
    for c in C:
        for h in H:
            prob += x_ch[c][h] == 1  
    
    return prob

prob= anadir_restriccion_xch(prob, x_ch, C, H)

def anadir_otras_restricciones(prob, H, C, e_ch, x_ch, E_c, beta, alpha, gamma_ch):
    ## Producción máxima
    for c in C:
        for h in H:
            prob += e_ch[c][h] <= x_ch[c][h] * E_c.loc[c][0]

    ## Producción mínima
    for c in C:
        for h in H:
            prob += e_ch[c][h] >= x_ch[c][h] * E_c.loc[c][0]*alpha[c]

    ## Cambio máximo de producción de energia
    for c in C:
        for h in range(1, len(H[0:-1])):
            prob += e_ch[c][H[h]] - e_ch[c][H[h-1]] <= beta[c] * E_c.loc[c][0]
            prob += e_ch[c][H[h-1]] - e_ch[c][H[h]] <= beta[c] * E_c.loc[c][0]    

    # Cambio de estado de cada central
    for c in C:
        for h in range(1, len(H[0:-1])):
            prob += gamma_ch[c][H[h]] == x_ch[c][H[h]] - x_ch[c][H[h-1]]  # Cambio de estado de encendido/apagado
            prob += gamma_ch[c][H[h]] == x_ch[c][H[h-1]] - x_ch[c][H[h]]  # Cambio de estado de encendido/apagado           

    return prob

prob= anadir_otras_restricciones(prob, H, C, e_ch, x_ch, E_c, beta, alpha, gamma_ch)

status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)
costo_inicial = value(prob.objective)

Status: Optimal

def almacenar_resultados_modelo(prob):   
    energy = {}
    x = {}
    gamma = {}
    result_plant = []
    result_energy = []
    result_hour = []
    block_size = 240
    block_count = 0
    current_dict = energy  # Start with energy

    for v in prob.variables():
        p = (v.name).split('_')[1]
        e = v.varValue
        h = (v.name).split('_')[2] + '-' + (v.name).split('_')[3] + '-' + (v.name).split('_')[4] + '-' + (v.name).split('_')[5]

        result_plant.append(p)
        result_energy.append(e)
        result_hour.append(h)

        block_count += 1
        if block_count % block_size == 0:
            if current_dict == energy:
                current_dict = x
            elif current_dict == x:
                current_dict = gamma
            else:
                break

            result_plant = []  
            result_energy = []
            result_hour = []

        current_dict['central'] = result_plant
        current_dict['hora'] = result_hour
        current_dict['energía'] = result_energy

    df_energy = pd.DataFrame(energy)
    df_xch = pd.DataFrame(x)
    df_gamma = pd.DataFrame(gamma)

    return (df_energy, df_xch, df_gamma)

df_result_energy, df_result_xch, df__result_gamma= almacenar_resultados_modelo(prob)

centrals = df_result_energy['central'].unique()

plt.figure(figsize=(20,10))

for central in centrals:
    plt.plot(df_result_energy[df_result_energy['central'] == central]['hora'], 
             df_result_energy[df_result_energy['central'] == central]['energía'], 
             marker='o', linestyle='-', 
             label=f'Energía producida por la central {central}')

    
plt.title(f'Energía producida por cada central por hora', weight= 'bold')
plt.xlabel('Hora')
plt.ylabel('Energía producida (MWh)')
plt.tick_params(axis='x', rotation=90)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()

centrals = df_result_energy['central'].unique()
fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    limite_superior = E_c.loc[central, 'producción']  # Obtiene el límite superior para la central
    limite_inferior = E_c.loc[central, 'producción'] * df_alpha.loc[central, 'alpha']  # Obtiene el límite inferior para la central
    axs[row, col].plot(df_result_energy[df_result_energy['central'] == central]['hora'], 
                       df_result_energy[df_result_energy['central'] == central]['energía'], 
                       marker='o', linestyle='-', 
                       label=f'Energía producida por la central {central}')
    axs[row, col].axhline(y=limite_superior, color='r', linestyle='--', label='Capacidad máxima de la central')
    axs[row, col].axhline(y=limite_inferior, color='b', linestyle='--', label='Mínima produccción admitida')
    axs[row, col].set_title(f'Producción de energia por central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Energía producida (MWh)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()

fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    df_central = df_result_energy[df_result_energy['central'] == central]
    df_central['cambio_energia'] = df_central['energía'].diff()  # Calcula la derivada de la energía
    
    limite_superior = E_c.loc[central, 'producción'] * df_beta.loc[central, 'beta']  # Obtiene el límite superior para la central
    limite_inferior = E_c.loc[central, 'producción'] * df_beta.loc[central, 'beta']  * -1 # Obtiene el límite inferior para la central
    
    axs[row, col].axhline(y=limite_superior, color='r', linestyle='--', label='Incremento de producción máximo admisible de la central')
    axs[row, col].axhline(y=limite_inferior, color='b', linestyle='--', label='Descenso de producción máximo admisible de la central')

    axs[row, col].plot(df_central['hora'], 
                       df_central['cambio_energia'], 
                       marker='o', linestyle='-', 
                       label=f'Cambio de energía en la central {central}')
    axs[row, col].set_title(f'Cambio de producción de energía por hora en la central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Cambio de energía producida (MWh/h)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()

df_result_energy['hora'] = pd.to_datetime(df_result_energy['hora'])
df_energia_demanda = pd.merge(df_result_energy, D_h, on='hora', how='left')
df_energia_demanda = df_energia_demanda.astype({'hora': 'string'})

produccion_total = df_energia_demanda.groupby(df_energia_demanda['hora'])[['energía', 'demanda']].sum().reset_index()
produccion_total['demanda'] = produccion_total['demanda']/10

plt.figure(figsize=(18,8))

plt.plot(produccion_total['hora'], produccion_total['energía'], label = 'Producción de energía total por hora', marker='x')

plt.xlabel('Hora')
plt.ylabel('MWh')
plt.xticks(rotation=90)

plt.legend()
plt.show

<function matplotlib.pyplot.show(close=None, block=None)>

plt.figure(figsize=(18,8))

plt.plot(produccion_total['hora'], produccion_total['energía'], label = 'Producción de energía total por hora', marker='x')
plt.plot(produccion_total['hora'], produccion_total['demanda'], label = 'Demanda de energía total por hora', linestyle = '--')

plt.fill_between(produccion_total['hora'], produccion_total['energía'], produccion_total['demanda'], 
                 where = (produccion_total['energía'] >= produccion_total['demanda']), facecolor='green', alpha = 0.5)

plt.xlabel('Hora')
plt.ylabel('MWh')
plt.xticks(rotation=90)

plt.legend()
plt.show

<function matplotlib.pyplot.show(close=None, block=None)>

# Calcular la diferencia entre la producción y la demanda
produccion_total['diferencia'] = produccion_total['energía'] - produccion_total['demanda']

plt.figure(figsize=(10,5))

plt.plot(produccion_total['hora'], produccion_total['diferencia'], label='Diferencia', color='purple')

plt.title('Diferencia entre la producción y la demanda')
plt.xlabel('Hora')
plt.ylabel('Diferencia (MWh)')
plt.legend()
plt.xticks(rotation=90)


plt.show()

# Transformamos el dataframe P_c a un formato largo
P_c_long = P_c.reset_index().melt(id_vars='index', var_name='hora', value_name='costo_combustible')
P_c_long.rename(columns={'index':'central'}, inplace=True)
P_c_long = P_c_long.astype({'hora': 'string'})

costos = pd.merge(df_energia_demanda, P_c_long,  how='left', left_on=['central','hora'], right_on = ['central','hora'])
costos['costo_combustible'] = costos['energía'] * costos ['costo_combustible']

O_c_long = O_c.reset_index().melt(id_vars='index', var_name='hora', value_name='costo_operativo')
O_c_long.rename(columns={'index':'central'}, inplace=True)
O_c_long = O_c_long.astype({'hora': 'string'})

costos = pd.merge(costos, O_c_long,  how='left', left_on=['central','hora'], right_on = ['central','hora'])

S_ch_T = S_ch.T
S_ch_T_long = S_ch_T.reset_index().melt(id_vars='central', var_name='hora', value_name='costo_sanitario')
S_ch_T_long.rename(columns={'index':'central'}, inplace=True)
S_ch_T_long = S_ch_T_long.astype({'hora': 'string'})

costos = pd.merge(costos, S_ch_T_long,  how='left', left_on=['central','hora'], right_on = ['central','hora'])
costos['costo_sanitario'] = costos['energía'] * costos ['costo_sanitario']

fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    df_central = costos[costos['central'] == central]

    df_central.groupby('hora').sum()[['costo_operativo', 'costo_combustible', 'costo_sanitario']].plot(kind='bar', stacked=True, ax=axs[row, col])
    
    axs[row, col].set_title(f'Distribución de costos de produccion de energía por hora en la central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Costo (Costo/MWh)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()

costos['costo total'] = costos['costo_operativo'] + costos['costo_combustible'] + costos['costo_sanitario']
costo_total = costos.groupby(costos['hora'])['costo total'].sum().reset_index()

# Crea un gráfico de barras con un ancho adecuado
plt.figure(figsize=(20, 6))
plt.bar(costo_total['hora'], costo_total['costo total'], label ='Costo total por hora')

plt.xlabel('Horas')
plt.ylabel('Costos')
plt.title('Costo Total por hora')
plt.legend()
plt.xticks(rotation=90)
plt.show()

costo_minimo = costo_total['costo total'].min()
costo_maximo = costo_total['costo total'].max()
hora_minima = costo_total['hora'].loc[costo_total['costo total']==costo_minimo].values[0]
hora_maxima = costo_total['hora'].loc[costo_total['costo total']==costo_maximo].values[0]

print(f'La hora menos costosa es: {hora_minima}')

La hora menos costosa es: 2023-07-01 09:00:00

print(f'La hora más costosa es: {hora_maxima}')

La hora más costosa es: 2023-07-01 20:00:00

costo_total_por_planta = costos.groupby(costos['central'])['costo total'].sum().reset_index()
plt.figure(figsize=(20, 6))
plt.bar(costo_total_por_planta['central'], costo_total_por_planta['costo total'], label ='Costo total por central')

plt.xlabel('Centrales')
plt.ylabel('Costos')
plt.title('Costo Total por central al día')
plt.legend()
plt.show()

costo_minimo = costo_total_por_planta['costo total'].min()
costo_maximo = costo_total_por_planta['costo total'].max()
central_min = costo_total_por_planta['central'].loc[costo_total_por_planta['costo total']==costo_minimo].values[0]
central_max = costo_total_por_planta['central'].loc[costo_total_por_planta['costo total']==costo_maximo].values[0]

print(f'La central menos costosa es: {central_min}')

La central menos costosa es: I

print(f'La central más costosa es: {central_max}')

La central más costosa es: A

coste_dia = value(prob.objective)
print(f'El coste total del día: {coste_dia}')

El coste total del día: 4508697.3013742585

def hallar_limite_viable(demanda, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha):
    prob = LpProblem('Problema_de_la_demanda', LpMinimize)
    e_ch, x_ch, gamma_ch = crear_variables_decision(C, H)
    prob = crear_funcion_obj(prob, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, C, H)
    prob = anadir_restriccion_xch(prob, x_ch, C, H)
    prob = anadir_otras_restricciones(prob, H, C, e_ch, x_ch, E_c, beta, alpha, gamma_ch)
    for h in H:
        prob+= (lpSum([e_ch[c][h] for c in C])) == demanda.loc[h][0]

    return prob, e_ch, x_ch, gamma_ch

coeficiente = 1
paso = 0.01
status = 'Optimal'

while (status == 'Optimal'):
    demanda = daily_demand * coeficiente
    prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(demanda, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)
    status = LpStatus[prob.solve()]
    coeficiente-= paso

porcentaje_inviable = 100 - (coeficiente+paso)*100
print(f"El porcentaje de reducción de demanda de energía que hace que la solución sea inviable es: ", porcentaje_inviable, "%")

El porcentaje de reducción de demanda de energía que hace que la solución sea inviable es:  58.00000000000005 %

limite_viable = daily_demand * (coeficiente + 2*paso)
prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(limite_viable, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)

status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)

Status: Optimal

910405.7955026062

df_result_limits, df_result_xch, df_result_gamma = almacenar_resultados_modelo(prob)
df_result_limits['hora'] = pd.to_datetime(df_result_limits['hora'])
df_energia_demanda_limite = pd.merge(df_result_limits, D_h, on='hora', how='left')
df_energia_demanda_limite = df_energia_demanda_limite.astype({'hora': 'string'})

produccion_total_limite = df_energia_demanda_limite.groupby(df_energia_demanda_limite['hora'])[['energía', 'demanda']].sum().reset_index()
produccion_total_limite['demanda'] = produccion_total_limite['demanda']/10

plt.figure(figsize=(18,8))

plt.plot(produccion_total_limite['hora'], produccion_total_limite['energía'], label = 'Producción de energía total por hora', marker='x')

plt.xlabel('Hora')
plt.ylabel('MWh')
plt.xticks(rotation=90)

plt.legend()
plt.show

<function matplotlib.pyplot.show(close=None, block=None)>

produccion_minima = produccion_total_limite['energía'].min()
hora_minima = produccion_total_limite['hora'].loc[produccion_total_limite['energía']==produccion_minima].values[0]
print(f'El momento en que la producción cae a su mínimo posible es: {hora_minima}')

El momento en que la producción cae a su mínimo posible es: 2023-07-01 08:00:00

demanda_minima = produccion_total_limite['demanda'].min()
hora_minima = produccion_total_limite['hora'].loc[produccion_total_limite['demanda']==demanda_minima].values[0]
print(f'El momento en que la demanda es menor es: {hora_minima}')

El momento en que la demanda es menor es: 2023-07-01 08:00:00

status = prob.solve()
print("Status: ", LpStatus[status])
costo_limite = value(prob.objective) 
value(prob.objective)

Status:  Optimal

910405.7955026062

dif =round((100 - ((costo_limite*100)/costo_inicial)),2)
print(dif)

79.81

print(f'El costo mínimo para este escenario es {costo_limite}, comparado con el costo inicial de {costo_inicial}, lo que indica una caída del {dif}% en el costo')

El costo mínimo para este escenario es 910405.7955026062, comparado con el costo inicial de 4508697.3013742585, lo que indica una caída del 79.81% en el costo

limite_viable = daily_demand * (coeficiente + paso)
prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(limite_viable, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)

status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)

Status: Infeasible

875612.3748042658

# Definicion del problema
prob = LpProblem('Optimizacion_Energia_ONOFF', LpMinimize)

# Debido a que el codigo del apartado anterior se definió modularmente y general, entonces:
#Creamos las variables decisión
e_ch, x_ch, gamma_ch = crear_variables_decision(C, H)

#No es necesario redifinir las variables, puesto que se mantiene el contexto.

# Definicion de la función objetivo:
prob= crear_funcion_obj(prob, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, C, H)

prob= anadir_restriccion_demanda(prob, e_ch, D_h, C, H) # Satisfacción de la demanda

def anadir_restriccion_xch2(prob, x_ch, C, H):
    centrales_nucleares = ['D', 'E']
    for c in C:
        for h in H:
            if (c in centrales_nucleares):
                prob += x_ch[c][h]==1 # Las centrales nucleares siempre deben estar encendidas
    return prob

prob = anadir_restriccion_xch2(prob, x_ch, C, H)

prob= anadir_otras_restricciones(prob, H, C, e_ch, x_ch, E_c, beta, alpha, gamma_ch)

status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)
costo_inicial_onoff = value(prob.objective)
costo_inicial_onoff

Status: Optimal

4497165.960429066

df_results, df_result_xch, df_result_gamma = almacenar_resultados_modelo(prob)

centrals = df_results['central'].unique()

plt.figure(figsize=(20,10))

for central in centrals:
    plt.plot(df_results[df_results['central'] == central]['hora'], 
             df_results[df_results['central'] == central]['energía'], 
             marker='o', linestyle='-', 
             label=f'Energía producida por la central {central}')

    
plt.title(f'Energía producida por cada central por hora', weight= 'bold')
plt.xlabel('Hora')
plt.ylabel('Energía producida (MWh)')
plt.tick_params(axis='x', rotation=90)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()

df_results = df_results.astype({'hora': 'datetime64[ns]'})
compared_results = pd.merge(df_results, df_result_energy, how= 'left', on=['central', 'hora'])

centrals = df_results['central'].unique()
fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    limite_superior = E_c.loc[central, 'producción']  # Obtiene el límite superior para la central
    limite_inferior = E_c.loc[central, 'producción'] * df_alpha.loc[central, 'alpha']  # Obtiene el límite inferior para la central
    axs[row, col].plot(df_results[df_results['central'] == central]['hora'], 
                       df_results[df_results['central'] == central]['energía'], 
                       marker='o', linestyle='-', 
                       label=f'Energía producida por la central {central}')
    axs[row, col].axhline(y=limite_superior, color='r', linestyle='--', label='Capacidad máxima de la central')
    axs[row, col].axhline(y=limite_inferior, color='b', linestyle='--', label='Mínima produccción admitida')
    axs[row, col].set_title(f'Producción de energia por central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Energía producida (MWh)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()

fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    df_central = compared_results[compared_results['central'] == central]

    axs[row, col].plot(df_central['hora'], 
                       df_central['energía_x'], 
                       marker='o', linestyle='-', 
                       label=f'Producción de energia con centrales encendidas y apagadas')
    
    axs[row, col].plot(df_central['hora'], 
                       df_central['energía_y'], 
                       marker='o', linestyle='--', 
                       label=f'Producción de energia con centrales siempre encendidas')
    
    axs[row, col].fill_between(df_central['hora'], df_central['energía_y'], df_central['energía_x'], 
                               where = (df_central['energía_y'] != df_central['energía_x']), facecolor= 'gray', alpha=0.5 )

    axs[row, col].set_title(f'Diferencia de energía producida en la central {central}, respecto al escenario anterior')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Energía producida (MWh/h)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()

centrals = df_result_xch['central'].unique()
fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2

    axs[row, col].plot(df_result_xch[df_result_xch['central'] == central]['hora'], 
                       df_result_xch[df_result_xch['central'] == central]['energía'], 
                       marker='o', linestyle='-', 
                       label=f'Energía producida por la central {central}')
    axs[row, col].axhline(y=1, color='r', linestyle='--', label='Capacidad máxima de la central')
    axs[row, col].axhline(y=0, color='b', linestyle='--', label='Mínima produccción admitida')
    axs[row, col].set_title(f'Producción de energia por central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Energía producida (MWh)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()

df_results = df_results.astype({'hora': 'string'})

new_costos = pd.merge(df_results, P_c_long, how='left', left_on=['central','hora'], right_on = ['central','hora'])
new_costos['costo_combustible'] = new_costos['energía'] * new_costos ['costo_combustible']

df_result_xch = df_result_xch.astype({'hora': 'datetime64[ns]'})
df_result_xch = df_result_xch.astype({'hora': 'string'})

new_costos = pd.merge(new_costos, O_c_long,  how='left', left_on=['central','hora'], right_on = ['central','hora'])
new_costos = pd.merge(new_costos, df_result_xch, how='left', on=['central', 'hora'])

new_costos['costo_operativo'] = new_costos['costo_operativo'] * new_costos['energía_y']

new_costos = pd.merge(new_costos, S_ch_T_long,  how='left', left_on=['central','hora'], right_on = ['central','hora'])
new_costos['costo_sanitario'] = new_costos['energía_x'] * new_costos ['costo_sanitario']

fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    df_central = new_costos[new_costos['central'] == central]

    df_central.groupby('hora').sum()[['costo_operativo', 'costo_combustible', 'costo_sanitario']].plot(kind='bar', stacked=True, ax=axs[row, col])
    
    axs[row, col].set_title(f'Distribución de costos de producción de energía por hora para la central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Costo (Costo/MWh)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()

new_costos['costo_total'] = new_costos['costo_operativo'] + new_costos['costo_combustible'] + new_costos['costo_sanitario']
new_costo_total = new_costos.groupby(new_costos['hora'])['costo_total'].sum().reset_index()

plt.figure(figsize=(20, 6))
plt.plot(costo_total['hora'], costo_total['costo total'], label ='Costo total por hora, anterior')
plt.plot(new_costo_total['hora'], new_costo_total['costo_total'], label ='Costo total por hora actual')

plt.xlabel('Horas')
plt.ylabel('Costos')
plt.title('Costo Total por hora')
plt.legend()
plt.xticks(rotation=90)
plt.show()

new_costo_total['diferencia'] = costo_total['costo total'] - new_costo_total['costo_total']

plt.figure(figsize=(10,5))

plt.plot(new_costo_total['hora'], new_costo_total['diferencia'], label='Diferencia', color='purple')

plt.title('Diferencia entre la producción y la demanda')
plt.xlabel('Hora')
plt.ylabel('Diferencia (MWh)')
plt.legend()
plt.xticks(rotation=90)


plt.show()

def hallar_limite_viable(demanda, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha):
    prob = LpProblem('Problema_de_la_demanda', LpMinimize)
    e_ch, x_ch, gamma_ch = crear_variables_decision(C, H)
    prob = crear_funcion_obj(prob, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, C, H)
    prob = anadir_restriccion_xch2(prob, x_ch, C, H)
    prob = anadir_otras_restricciones(prob, H, C, e_ch, x_ch, E_c, beta, alpha, gamma_ch)
    for h in H:
        prob+= (lpSum([e_ch[c][h] for c in C])) == demanda.loc[h][0]

    return prob, e_ch, x_ch, gamma_ch

coeficiente = 1
paso = 0.005
status = 'Optimal'

while (status == 'Optimal'):
    demanda = daily_demand * coeficiente
    prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(demanda, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)
    status = LpStatus[prob.solve()]
    coeficiente+= paso

porcentaje_inviable = ((coeficiente+paso)*100)-100
print(f"El porcentaje de reducción de demanda de energía que hace que la solución sea inviable es: ", porcentaje_inviable, "%")

El porcentaje de reducción de demanda de energía que hace que la solución sea inviable es:  32.99999999999929 %

limite_viable = daily_demand * 1.315
prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(limite_viable, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)

status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)

Status: Optimal

7755835.807726487

df_result_limits, df_result_xch, df_result_gamma = almacenar_resultados_modelo(prob)
df_result_limits['hora'] = pd.to_datetime(df_result_limits['hora'])
df_energia_demanda_limite = pd.merge(df_result_limits, D_h, on='hora', how='left')
df_energia_demanda_limite = df_energia_demanda_limite.astype({'hora': 'string'})

df_result_limits, df_result_xch, df_result_gamma = almacenar_resultados_modelo(prob)
df_result_limits['hora'] = pd.to_datetime(df_result_limits['hora'])
df_energia_demanda_limite = pd.merge(df_result_limits, D_h, on='hora', how='left')
df_energia_demanda_limite = df_energia_demanda_limite.astype({'hora': 'string'})

produccion_total_limite = df_energia_demanda_limite.groupby(df_energia_demanda_limite['hora'])[['energía', 'demanda']].sum().reset_index()
produccion_total_limite['demanda'] = produccion_total_limite['demanda']/10

plt.figure(figsize=(18,8))

plt.plot(produccion_total_limite['hora'], produccion_total_limite['energía'], label = 'Producción de energía total por hora', marker='x')

plt.xlabel('Hora')
plt.ylabel('MWh')
plt.xticks(rotation=90)

plt.legend()

<matplotlib.legend.Legend at 0x213fd514a40>

produccion_minima = produccion_total_limite['energía'].max()
hora_minima = produccion_total_limite['hora'].loc[produccion_total_limite['energía']==produccion_minima].values[0]
print(f'El momento en que la producción llega a su máximo posible es: {hora_minima}')

El momento en que la producción llega a su máximo posible es: 2023-07-01 20:00:00

demanda_minima = produccion_total_limite['demanda'].max()
hora_minima = produccion_total_limite['hora'].loc[produccion_total_limite['demanda']==demanda_minima].values[0]
print(f'El momento en que la demanda es mayor es: {hora_minima}')

El momento en que la demanda es mayor es: 2023-07-01 20:00:00

status = prob.solve()
print("Status: ", LpStatus[status])
costo_limite = value(prob.objective) 
value(prob.objective)

Status:  Optimal

7755835.807726487

dif =round((((costo_limite*100)/costo_inicial)),2)-100
print(dif)

72.02000000000001

print(f'El costo máximo para este escenario es {costo_limite}, comparado con el costo inicial de {costo_inicial_onoff}, lo que indica una subida del {dif}% en el costo')

El costo máximo para este escenario es 7755835.807726487, comparado con el costo inicial de 4497165.960429066, lo que indica una subida del 72.02000000000001% en el costo

limite_viable = daily_demand * (coeficiente + paso)
prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(limite_viable, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)

status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)

Status: Infeasible

7941831.5364006655

	demanda
hora
2023-07-01 01:00:00	22424.578280
2023-07-01 02:00:00	20288.860640
2023-07-01 03:00:00	17962.060638
2023-07-01 04:00:00	16585.284168
2023-07-01 05:00:00	15976.342990
...	...
2023-07-31 20:00:00	27991.637100
2023-07-31 21:00:00	27942.225340
2023-07-31 22:00:00	27630.460640
2023-07-31 23:00:00	27299.872400
2023-08-01 00:00:00	26400.119460

	tipo	combustible	producción
central
A	carbón	carbón sub-bituminoso	5720.22000
B	carbón	carbón bituminoso	5244.74856
C	carbón	carbón bituminoso	861.36600
D	nuclear	nuclear	3252.82456
E	nuclear	nuclear	4177.39200
F	petróleo	licor negro	151.67496
G	biomasa	licor negro	941.91944
H	gas	gas natural	3522.61048
I	hidráulica	agua	536.11236
J	carbón	carbón bituminoso	8793.66424

	coste_comb	coste_oper	coste_on_off
combustible
agua	0.000	8.176099	0.0
carbón bituminoso	18.050	7.384495	147.0
carbón sub-bituminoso	9.195	7.380880	147.0
gas natural	17.435	5.287164	72.0
licor negro	7.440	5.295063	72.0
nuclear	3.615	16.310160	0.0

		coste_salud
hora	central
2023-07-01 01:00:00	A	2.148816
	B	9.030773
	C	3.422946
2023-07-01 02:00:00	A	2.022208
2023-07-01 02:00:00	B	7.706546
...	...	...
2023-07-31 23:00:00	B	13.274831
2023-07-31 23:00:00	C	14.214563
2023-08-01 00:00:00	A	2.562029
	B	12.606380
	C	10.906997

Optimizando la generación de electricidad¶

Datos¶

Estimación de demanda¶

Centrales¶

Costes¶

Restricciones y suposiciones¶

Optimización¶

Apartado A (4 puntos)¶

1. Formulación matemática¶

1.1 Esquema visual:¶

Definición de variables:¶

2. Formulación en la librería PulP¶

Definición de variables¶

Definición de la función objetivo¶

Definición de las restricciones¶

3. Resolver el problema y analizar la solución obtenida, pintando:¶

Cuánto produce cada central en cada hora.¶

¿Hay alguna que produzca lo mínimo posible? Si es así, ¿durante qué horas del día?, ¿Hay alguna que produzca lo máximo posible? Si es así, ¿durante qué horas del día?¶

¿Hay alguna que cambie entre una hora y la siguiente lo máximo posible? Si es así, ¿creciendo o decreciendo de producción?¶

Cuál es la producción total por hora¶

¿Cómo compara con respecto a la demanda?¶

Cuál es el coste incurrido por cada central en cada hora, desglosado en sus componentes (combustible, operación y sanitario).¶

Cuál es el coste total por hora.¶

¿Qué horas son las más y las menos costosas?¶

Cuál es el coste por central en el día.¶

¿Qué centrales son las más y las menos costosas?¶

Cuál es el coste total del día¶

4. ¿Cuánto tiene que caer la demanda en % para que el problema esté en el límite de ser o no factible? Aquí puedes suponer que la demanda decae en el mismo porcentaje para las 24 horas del día. Muestra que para ese valor límite:¶

En algún momento del día la producción cae a su mínimo posible. ¿Qué momento es ese y por qué?¶

El coste es más bajo que en el escenario anterior. ¿Cuál sería este coste mínimo?¶

Si decae todavía un poco más, el problema se convierte en no factible.¶

Apartado B (6 puntos)¶

1. Formular el problema matemáticamente de cómo se debe operar a lo largo del dia¶

2. Trasladar esta nueva formulación a la librería PuLP.¶

Variables de decisión¶

Función Objetivo¶

Definición de las restricciones¶

3. Resolver el problema y analizar la solución obtenida, comparando con el resultado obtenido antes:¶

¿Cómo varía la producción de energía de cada central? La demanda sigue siendo la misma, pero es posible que no se opere exactamente igual que antes..¶

¿Hay alguna central que se desactive o active? Si es así, ¿en qué horas?¶

¿Hay alguna siempre encendida (además de las nucleares)?¶

¿Y alguna siempre apagada?¶

¿Cómo varía la distribución de costes? ¿Es más barato o más caro que antes el total del día? ¿Por qué?¶

4. ¿Cuánto tiene que subir la demanda en % para que el problema esté en el límite de ser o no factible? De nuevo supón que la demanda aumenta en el mismo porcentaje para las 24 horas del día. Muestra que para ese valor límite:¶

Apartado C (opcional)¶

Criterios de evaluación¶

2. Formulación en la librería `PulP`¶

2. Trasladar esta nueva formulación a la librería `PuLP`.¶